#ifndef _FSES3_H_
#define _FSES3_H_

#include "Equations/TwoPhaseEquation.hpp"

namespace Tuna {

  template <typename T, int Dim> class FSES3;

  template <typename T>
  struct Typeinfo <FSES3<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };

  template <typename T>
  struct Typeinfo <FSES3<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };

  template <typename T>
  struct Typeinfo <FSES3<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };
  
  template<typename Tprec, int Dim>
  class FSES3 : public TwoPhaseEquation<FSES3<Tprec, Dim> >
  {

    typedef TwoPhaseEquation<FSES3<Tprec, Dim> > TP_FSES3;

    using GeneralEquation< TP_FSES3 >::aE;
    using GeneralEquation< TP_FSES3 >::aW;
    using GeneralEquation< TP_FSES3 >::aN;
    using GeneralEquation< TP_FSES3 >::aS;
    using GeneralEquation< TP_FSES3 >::aF;
    using GeneralEquation< TP_FSES3 >::aB;
    using GeneralEquation< TP_FSES3 >::aP;
    using GeneralEquation< TP_FSES3 >::sp;
    using GeneralEquation< TP_FSES3 >::dx;
    using GeneralEquation< TP_FSES3 >::dy;
    using GeneralEquation< TP_FSES3 >::dz;
    using GeneralEquation< TP_FSES3 >::dt;
    using GeneralEquation< TP_FSES3 >::bi;
    using GeneralEquation< TP_FSES3 >::ei;
    using GeneralEquation< TP_FSES3 >::bj;
    using GeneralEquation< TP_FSES3 >::ej;
    using GeneralEquation< TP_FSES3 >::bk;
    using GeneralEquation< TP_FSES3 >::ek;
    using GeneralEquation< TP_FSES3 >::applyBoundaryConditions1D;
    using GeneralEquation< TP_FSES3 >::applyBoundaryConditions2D;
    using GeneralEquation< TP_FSES3 >::applyBoundaryConditions3D;

    using TP_FSES3::S;
    using TP_FSES3::p;
    using TP_FSES3::phi_0;
    using TP_FSES3::Srw;
    using TP_FSES3::Sro;
    using TP_FSES3::mu_w;
    using TP_FSES3::mu_o;
    using TP_FSES3::k;
    using TP_FSES3::porosity;
    using TP_FSES3::injection;

  public:
    typedef Tprec prec_t;
    typedef typename TunaArray<prec_t, Dim >::huge ScalarField;
    
    FSES3() : TwoPhaseEquation<FSES3<prec_t, Dim > >() { }    
    ~FSES3() { };
    
    inline bool calcCoefficients1D(); 
    inline bool calcCoefficients2D() { };
    inline bool calcCoefficients3D();
    inline void printInfo() { std::cout << " FSES3 "; }
  };

/*
 *  Lineal for realtive permeability , Average for Sw
 */
template<typename Tprec, int Dim>
inline bool FSES3<Tprec, Dim>::calcCoefficients1D () 
{
    static prec_t Sw_e, Sw_w;

    // Lineal 
    static prec_t mult = k * dt / 
      ( porosity * dx * dx * (1 - Srw - Sro) * mu_w );
      
    aE = 0.0; aW = 0.0; aP = 0.0; sp = 0.0;

    for (int i =  bi; i <= ei; ++i) {      

      // Average
      Sw_e = ( phi_0(i+1) + phi_0(i) ) * 0.5;
      Sw_w = ( phi_0(i-1) + phi_0(i) ) * 0.5; 

      // Lineal
      aE (i) = (Sw_e - Srw) * mult;
      aW (i) = (Sw_w - Srw) * mult;
      aP (i) = aE (i) + aW (i);      
    }

    // Dirichlet right side
    aP(ei) += aE(ei);
    sp(ei) = 2 * aE(ei) * p(ei+1);
    aE(ei) = 0;

    // Neumann left side
    aP(bi) -= aW(bi) ;
    sp(bi) = aW(bi) * dx * ( 3.47e-7 * mu_w / k) ;
    aW(bi) = 0;

    return 0;
}


/*
 *  Lineal for realtive permeability , Upwind for Sw
 */
template<typename Tprec, int Dim>
inline bool FSES3<Tprec, Dim>::calcCoefficients3D () 
{
  static prec_t Sw_e, Sw_w, Sw_n, Sw_s, Sw_f, Sw_b;

    // Lineal 
  static prec_t multx = k * dt / 
    ( porosity * dx * dx * (1 - Srw - Sro) * mu_w );

  static prec_t multy = k * dt / 
    ( porosity * dy * dy * (1 - Srw - Sro) * mu_w );
 
  static prec_t multz = k * dt / 
    ( porosity * dz * dz * (1 - Srw - Sro) * mu_w );

  aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aF = 0.0; aB = 0.0; aP = 0.0; sp = 0.0;

  for (int ki = bk; ki <= ek; ++ki) 
    for (int i =  bi; i <= ei; ++i)
      for (int j = bj; j <= ej; ++j) {  

      // Average
	Sw_e = ( phi_0(i+1,j,ki) + phi_0(i,j,ki) ) * 0.5;
	Sw_w = ( phi_0(i-1,j,ki) + phi_0(i,j,ki) ) * 0.5;
	Sw_n = ( phi_0(i,j+1,ki) + phi_0(i,j,ki) ) * 0.5;
	Sw_s = ( phi_0(i,j-1,ki) + phi_0(i,j,ki) ) * 0.5;
	Sw_f = ( phi_0(i,j,ki+1) + phi_0(i,j,ki) ) * 0.5;
	Sw_b = ( phi_0(i,j,ki-1) + phi_0(i,j,ki) ) * 0.5;

      // Lineal
	aE (i, j, ki) = (Sw_e - Srw) * multx;
	aW (i, j, ki) = (Sw_w - Srw) * multx;
	aN (i, j, ki) = (Sw_n - Srw) * multy;
	aS (i, j, ki) = (Sw_s - Srw) * multy;
	aF (i, j, ki) = (Sw_f - Srw) * multz;
	aB (i, j, ki) = (Sw_b - Srw) * multz;
	aP (i, j, ki) = aE (i, j, ki) + aW (i, j, ki) + 
	  aN (i, j, ki) + aS (i, j, ki) + aF (i, j, ki) + aB (i, j, ki);
    }    

  /* ----- Boundary conditions ----- */

  applyBoundaryConditions3D();


  // Source for injection : 
  sp (bi, bj, bk) +=  injection * dt / (porosity * dx * dy * dz);

/*
  sp (ei, ej) +=  injection * dt / (porosity * dx * dy);
  sp (bi, ej) +=  injection * dt / (porosity * dx * dy);
  sp (ei, bj) +=  injection * dt / (porosity * dx * dy);
*/

  // Sink for extraction :  
  sp (ei, ej, ek) -= injection * dt * phi_0(ei, ej, ek) / ( porosity * dx * dy * dz);
  //  sp (ei, ej, ek) -= injection * dt / ( porosity * dx * dy * dz);

   /***
    // I'm using an explicit scheme, so the boundary conditions are
    // different than in the implicit case.

  Range I(bi, ei), J(bj, ej), K(bk, ek);

    // Dirichlet right side
  aP(ei, J, K) += aE(ei, J, K);
  sp(ei, J, K) = 2 * aE(ei, J, K) * p(ei+1, J, K);
  aE(ei, J, K) = 0;

  // Neumann left side
  // Remember that Sw_e = 0.8 on the boundary, that is
  // the reason why sp(bi) has the form below.
  aP(bi, J, K) -= aW(bi, J, K) ;
  sp(bi, J, K) = aW(bi, J, K) * dx * ( 3.47e-7 * mu_w / k) ;
  aW(bi, J, K) = 0;
   ***/
    return 0;
}



} // Tuna namespace


#endif //_FSES3_H_

















